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ABSTRACT 

Recent pulsar surveys have increased the number of observed double neutron stars (DNS) in our 
galaxy enough so that observable trends in their properties are starting to emerge. In particular, 
it has been noted that the majority of DNS have eccentricities less than 0.3, which are surprisingly 
low for binaries that survive a supernova explosion that we believe imparts a significant kick to the 
neutron star. To investigate this trend, we generate many different theoretical distributions of DNS 
eccentricities using Monte Carlo population synthesis methods. We determine which eccentricity 
distributions are most consistent with the observed sample of DNS binaries. In agreement with 
Chaurasia & Bailes (2005), assuming all double neutron stars are equally as probable to be discovered 
as binary pulsars, we find that highly eccentric, coalescing DNS are less likely to be observed because 
of their accelerated orbital evolution due to gravitational wave emission and possible early mergers. 
Based on our results for coalescing DNS, we also find that models with vanishingly or moderately small 
kicks (a < 50 kins" 1 ) are inconsistent with the current observed sample of such DNS. We discuss the 
implications of our conclusions for DNS merger rate estimates of interest to ground-based gravitational- 
wave interferometers. We find that, although orbital evolution due to gravitational radiation affects 
the eccentricity distribution of the observed sample, the associated upwards correction factor to merger 
rate estimates is rather small (typically 10-40%). 

Subject headings: binaries: general — binaries: close — stars: neutron — pulsars: general — gravita- 
tional waves 



1. INTRODUCTION 

Since the discovery of PSR B1913 +16, the Hulse 
fc T aylor binary pulsar in 1974 (Hul se fc Tavlorl 
|1975[K seven more double neutron stars (DNS) 
have been found. Of all these, five sys- 

tems are considered "close" (or "coal es cing") : 
PSRs B1913 +16, B1534+12 iWolszczanl [l991), 
B212 7+11C llPhinnev fc Sigurdssonl I1991D. J1756- 
2251 l|Faulkner et al]l2004lL .I0737-3039A (|Burgav et all 
|2003|) . Three are "wide" (wi ll not coalesce in a Hubble 
time): PSRs J1518+4 904 iNice et all ll996T). .11811- 
1736 l|Lvne et all 12000). .11829+2456 l|Chamnion'et~aTI 
2004). A more recently disco vered binary pulsar 
J1906+0746 ijLorimer et alJ 12006ft is another possible 
close DNS, although its nature as such is not yet 
confirmed. Binaries are characterized as either close 
or wide depending on their evolutionary timescales. A 
DNS goes through an inspiral phase and eventually 
merges because energy is taken away from its orbit in 
the form of gravitational radiation. If the time between 
the birth of the primordial binary and the merger of the 
resulting DNS is less than a Hubble time, we consider it 
to be "coalescing" (close), and wide if not. Coalescing 
DNS are of interest to LIGO and other ground-based 
interferometers because DNS mergers are primary 
targets for gravitational-wave detectors. The relativistic 
nature of such inspiraling DNS make them objects of 
extreme interest for astrophysics and relativity studies. 



TABLE 1 
Observed DNS a Binaries 



PSR 


e 


a (R Q ) 


P orb (days) 


B1913+16 


0.617 


1.009 


0.323 


B1534+12 


0.274 


1.606 


0.421 


J1756-2251 


0.181 


1.187 


0.320 


J0737-3039 


0.088 


0.610 


0.102 


J1906+0746 


0.085 


0.612 


0.166 


J1518+4904 


0.249 


8.634 


8.634 


J1811-1736 


0.828 


14.982 


18.779 


J1829+2456 


0.139 


3.117 


1.176 
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Tombaugh Fellow 



a Eccentricity, projected semimajor axis, and orbital period of ob- 
served DNS (and recently discovered binary J1906+0746). Infor- 
mation obtained from the Australia Telescope National Facility 
(ATNF) Pulsar Catalogue. 



Of the coalescing DNS, PSR B2127-11C is a system 
with an e ccentricity of 0.68, located in the globular clus- 
ter M15 l|Phinnev fc Sigurdssonl IT99l|) . Star formation 
in the globular cluster stopped long before the birth 
of the recycled pulsar. This particular binary was al- 
most certainly created as the result of cluster dymaics. 
Our calculations in the present paper do not account 
for DNS formation in clusters, and therefore we will 
not include PSR B2127-11C in our analysis. Although 
PSR J1906+0746 is not yet a confirmed DNS, we include 
it in our analysis as such. 

The orbital properties of the binaries we consider here 
are given in Table ^ It is evident that six of the 
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eight observed systems have orbital eccentricities below 
0.3. The noticeable trend towards low eccentricity DNS 
is puzzling; one would think that a binary surviving 
two supernovae, belie ved to impart velo city kicks of up 
to lOOOkms" 1 (e.g., iHobbs et all 120051 and references 
therein), would typically be highly eccentric. Part of 
this can be explained by circularization that occurs af- 
ter the binary's first supernova. Also, a class of low ec- 
centricity DNS, which receive low velocity kicks at the 
bir th of the second ne utron star, has been proposed 
by Ivan den Heuvell l)2004f) , as a possible explanation for 
the low eccentrici ties (see also lPo^ siadlowski et al1l2004i: 
iPfahl et alJl2002l iDewi et al.ll2005lL However, evidence 
for moderate to large kicks exists when when we consider 
the distribution of isolated pulsar velocities ( Hobbs et al. 
120051 lArzoumanian et al.ll2002t iLvne fc Lorimer| ^994'l . 

In this paper we examine theoretical distributions 
of DNS eccentricities based on various model assump- 
tions. These distributions are generated by Monte 
Carlo methods using binary evolution synthesis mod- 
els (§2). We analyze the resulting eccentricity distri- 
butions from a number of different models (§3), paying 
special attention to small vs. standard kicks and birth 
vs. present populations. We perform Bayesian statis- 
tical tests to evaluate whether any and which models 
are more likely given the observed DNS sample (§4). 
In § 5, we show that high-eccentricity DNS are asso- 
ciated with shorter merger times. Systems with short 
merger times are harder to detect because of their early 
merg ers (see also lGondek-Rosinska. Bulik. fc Belczvnskl 
2005). Therefore high-eccentricity DNS are not likely 
to contribute to the observed sample, as also concluded 
bv lChaurasia fc Bailesl i|2005fi 2 . We conclude with a dis- 
cussion of our findings and their implications for DNS 
inspiral detections by ground-based gravitational-wave 
interferometers. 

2. THEORETICAL METHODS 
2.1. StarTrack: Monte Carlo Population Synthesis 

To understand the observed sample of DNS and their 
eccentricities, we must provide a context for these ob- 
servations by understanding the Galactic DNS popula- 
tion and their properties. Theoretical models of the 
stellar population are generated by population synthe- 
sis methods. Our population synthesis utilizes Monte 
Carlo techniques that randomly choose stellar parame- 
ters drawn from assumed distribution functions of ini- 
tial conditions. These initial parameters characterize a 
primordial binary, and from these initial conditions the 
evolution of the binary is followed, based on best cur- 
rent our understanding of astrophysical processes. This 
is done for a large number of binaries, enough so that an 
adequate number of objects of interest form (in this case, 
DNS) for us to examine statistically. 

We use the well-tested StarTrack Population Synthesis 
Code to generate a large number of binaries that sim- 
ulates the population and evolution of binaries in our 
galaxy. A more complete descri ption o f the code can be 
found in lBelczvnski et all (|2002fl and in lBelczvnski et all 
(2005). Monte Carlo techniques are used to draw the 

2 We note that lUh aurasia & Baile^ 120051) appeared on the 
LANL preprint server when we were already preparing this 
manuscript for submission to the Journal. 



initial parameters of the primordial binaries (masses, or- 
bital separation, and eccentricity) from specified distri- 
butions, which are described in the next section. Further- 
more the following are chosen from appropriate distribu- 
tions: (a) the time of birth for each primordial binary, 
chosen from 0-15 Gyr, with a constant star formation 
rate assumed, (b) the direction and magnitude of the 
kick imparted to a neutron star (NS) or black hole (BH) 
in an asymmetric supernova (SN) explosion, and (c) the 
position in the orbit where the SN takes place (relevant 
to eccentric orbits before the explosion). Other aspects 
of the binary's evolution are based on physical laws that 
depend on the binary's properties with some necessary 
parametrizations. 

We track the evolution of each binary and its physi- 
cal properties at appropriate time steps from the zero- 
age main sequence until the formation of compact ob- 
jects (white dwarfs, neutron stars, and black holes) and 
their subsequent orbital evolution. Effects such as stellar 
winds, stable and unstable mass transfer through Roche- 
lobe overflow and winds, angular momentum losses, su- 
pernovae, tidal interactions of binary components, and 
orbital decay due to the emission of gravitational radia- 
tion are taken into account. The evolution of a binary 
ends when the binary merges, or when the end time of 
the simulation is reached (in this case, at 15 Gyr). 

2.2. Model Assumptions 

For this paper, we look at ten different population syn- 
thesis models, denoted as A, C, El, Gl, G3, H2, LI, 
L2, Ml, and M2 [the naming convention originates in 
iBelczvnski et alJ (|2002)]. Model A is our reference or 
standard model described next. All the other models are 
derived from this, with one of the synthesis model pa- 
rameters varied. 

The standard model draws masses for the primordial 
primary (the more massive binary component) from the 
initial mass function ^>(M) oc M~ 2 - 7 in the mass range 
from 4Af Q to 100M Q , the masses relevant for neutron 
star and black hole formation. The binary mass ratio 
q, defined as the ratio of the smaller mass to the larger 
such that < q < 1, is drawn from a flat distribution, 
$(g) = 1, so that given a primary mass Mi, the mass of 
the secondary, M 2 , is M2 = qM\. The initial separation 
of a binary is chosen from a distribution that is assumed 
to be logarithmically flat: T(A) oc -j. Initial eccentrici- 
ties follow the thermal-equilibrium eccentricity distribu- 
tion, S(e) = 2e. For cases of dynamically unstable mass 
transfer and common-envelope evolution, we adopt an 
effective efficiency of otcE x A = 1; neutron stars enter- 
ing a common envelope phase are modeled taking into 
account hypercritical accretion. Helium star evolution 
is mod eled in accordance to the results bv llvanova et alJ 
(2003). Mass transfer phases between non-degenerate 
stars are assumed to be non-conservative (50% of the 
transfered mass is lost from the binary) . This lost mass is 
assumed to carry a specific angular momentum equal to 
2irjA 2 /P (i.e., j — 1). For supernova kicks imparted to 
neutron stars at birth we ado pt the joint Maxwell i an ve- 
locity distribution derived bv lArzoumanian et al.1 (2002) 
with characteristic velocities a\ = 90 km/s (weighted at 
40%) and cr 2 = 500 km/s (weighted at 60%). Hereafter 
we refer to this distribution of kicks as the "Arzoumanian 
kick distribution" . 
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TABLE 2 
Model Assumptions 



Model Description 

A Standard (see text for details) 

C No hypercritical accretion onto NS/BH 

in common envelopes 

El Effective common-envelope efficiency 
acE X A = 0.1 

G1,G3 Stellar winds reduced by f w i„d = 0.5,0.2 

H2 Helium giants do not develop 
deep convective envelopes 

L1,L2 Angular momentum of material lost in 

non-conservative mass transfer: j = 0.5, 2.0 

M1,M2 Initial mass ratio distribution: &(q) <x q~ 2 7 ,q 3 



The distinguishing characteristics of the other 9 mod- 
els are described in Table [21 With this set of models 
we vary the values (within reasonable ranges) of a num- 
ber of stellar and binary evolution parameters that are 
known to affect DNS formation significantly. For more 
information on the specific choi ces of values for both th e 
standard and these models, see IBelczvnski elaH {2002). 

As a possible explanation of the low eccentricities, it 
has been suggested that a class of neutron stars form- 
ing in DNS systems that received small-magnitudes kicks 
is responsible for th e obser ved low eccentricities (see 
van den Heuve]|l2004 see also IPodsiadlowski et alJl2004t 
Pfahl et alJl200a iDewi et al.ll2005|) . To investigate this 
suggestion statistically we created another 3 sets of mod- 
els which each adopt a different NS kick distribution 
(other than the Arzoumanian kick distribution): (i) a 
set of (unphysical in our opinion) models with zero NS 
kicks (all SN events are assumed to be symmetric); (ii) a 
set with a single Maxwellian of a = lOkrns -1 ; and (hi) 
a set with a single Maxwellian of a — 50kms _1 . With 
these choices we hope to examine the effect of zero or 
low kick magnitudes on DNS orbital eccentricities. Each 
of these three kick distributions are adopted for all ten 
models listed in Table [3 

We note that the sugges t ion in favor of low kicks ma de 
by l)van den Heuvell 120041: IPodsiadlowski et alJl200l|) is 
connected to NS forming from helium-rich progenitors 
with masses up to ~ 3.5 M Q (H-ric h progenitors in the 
mass range 8—13 M Q ; for details see lPodsiadlowski et alJ 
|2004|) . Examination of the helium-rich DNS progenitors 
in all our simulations reveals that the vast majority (typ- 
ically above 90% and in some cases close to 100%) have 
masses < 3.5 M©. Therefore our consideration of models 
with lo w kicks for all N S is co nsist ent with the sugges- 
tions b v |van den Heuvell l|2004fl and IPodsi adlowsk i et alJ 
(2004). Furthermore the models with perfectly symmet- 
ric explosions provide a clean case that allows us to track 
the qualitative behavior of DNS eccentricities to pre-SN 
properties, as discussed in what follows. 

3. ECCENTRICITY DISTRIBUTIONS 
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Fig. 1. — Eccentricity distribution functions for coalescing bina- 
ries at birth for three of the ten models in Table and for three 
different NS kick assumptions: Arzoumanian kicks (solid lines), 
Maxwellian with a = 50kms _1 (dashed lines), and zero kicks (dot- 
ted lines). Only close binaries are considered. 



From the large populations of binaries we evolve, we 
consider all of the DNS systems formed and focus on 
the resultant DNS eccentricity distrbution functions, for 
each model. We typically generate ~ 30 DNS for every 
10 5 primordial binaries evolved, though in some cases 10 5 
primordials will only produce a few relevant objects. In 
short, generating reliable statistics for a variety of dif- 
ferent models is computationally expensive and requires 
much time and patience. We make sure that for each 
model considered here we typically have model samples 
of 2, 000 - 3, 000 DNS at birth and at least ~ 200 DNS 
for models that have atypically low DNS formation rates 
(e.g., L2, Ml). 

3.1. Zero, Low, and "Standard" Neutron- Star Kicks: 
Eccentricity Distributions at Birth 

In Figs.n an d0we present the eccentricity distribution 
functions at DNS birth for three of the models considered 
in our study, and for close and wide DNS, respectively. 
Before we proceed with our statistical analysis and com- 
parison to the observed sample, there are a number of 
qualitative conclusions we can draw from these distribu- 
tions. 

In the case of coalescing DNS (Fig-QJ, we find that the 
eccentricity distributions are distinctly different for the 
models with the Arzoumanian ("standard") kick distri- 
bution and those with zero or low kicks. In the first case 
the distribution covers the full range of values up to unity, 
and for most models it is almost flat across the range. As 
kick magnitudes decrease the distributions appear to be 
restricted to a low range of values, most typically below 
e = 0.6. This behavior appears to be consistent across 
all ten binary evolution models, independent of the spe- 
cific model assumptions, an d it is qualitatively cons istent 
with the suggestion made bv lvan den Heuvell ()2004j) that 
very small NS kicks could be responsible for the low ec- 
centricities in the observed sample. 

When supernova explosions are perfectly symmetric or 
have very small kicks (~ 10 — 20kms~ 1 ) DNS eccen- 
tricities are determined by mass loss at collapse. The 
qualitative behavior seen in the results of our numerical 
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eccentricity 

Fig. 2. — Eccentricity distribution functions for wide binaries at 
birth. Line types as in Fig. 1. 

simulations can be understood based on the analytical 
relations that govern the orbital dynamics in this case. 
Following Kal ogeral l|1996|h we characterize the mass loss 
at NS formation using the parameter (3 defined as: 

_ Mns + M 2 
P ~ M 1 + M 2 ' [ ' 

where Mi and M 2 are the masses of the exploding pro- 
genitor and the companion, respectively, just prior to the 
explosion, and M/vs is the mass of the resulting neutron 
star. 

Then the eccentricity after symmetric explosions as- 
suming circular pre-SN orbits is 



It is in the models of purely symmetric or very small kick 
magnitudes that the distribution of helium-star progen- 
itors to coalescing DNS is revealed. 

The result of eccentricities being constrained to val- 
ues smaller than 0.6 clearly implies tightly constrained 
amounts of supernova mass loss as well, for these mod- 
els: (3 > 0.6, i.e., up to 40% mass loss relative to the 
pre-SN total binary mass. Such a constraint is con- 
firmed when we examine our numerical results for /?. 
For NS masses that are all equal to a constant value 
(e.g., I.4M0), this constraint would imply a very nar- 
row range of helium-mass progenitors (< 2.2 Mq). Such 
a tight upper limit on NS progenitor masses is marginally 
consistent with stellar evolution models that in dicate a 
minim um of 2.1 M Q for electron-capture collapse ijHabetsI 
1986). For the wider range of NS masses considered in 
our simulations (up to 2M Q ), we find that helium-star 
progenitors up to 3-4 M become possible. This mass 
distribution is determined by the binary evolution prior 
to the supernova event leading to DNS formation; as such 
we find that it is not greatly affected by the presence of 
significant kicks or not. Their presence though affects the 
DNS eccentricities; the stronger kicks are relative to the 
pre-SN orbital velocities, the more eccentricities spread 
to values over the whole range from to 1 and the one- 
to-one connection between eccentricity and SN mass loss 
is lost. 

For wide DNS (Fig.[2J), the behavior of the eccentricity 
distributions is much less dependent on the presence of 



kicks. We find that the SN mass loss amounts cover a 
wider range with helium-star progenitors as massive as 
8 — 10 M Q . As a result, eccentricities cover the full range 
up to unity with no apparent qualitative difference or 
bias towards low values. 

It is also interesting to note that for both close and 
wide binaries, in the case of zero NS kicks a double-peak 
structure is most clearly evident. Further investigation 
of this structure reveals that it is due to a double-peak 
structure in the iron core masses of helium stars, the 
NS progenitors (as first revealed and d iscussed in detail 
bv lTimmes. Wooslev. fc Weaverl li"996|) . This structure 
propagates to the eccentricity through the mass-loss dur- 
ing the second supernova explosion in the systems: the 
range of helium-star masses is narrow, and therefore the 
double-peak structure persists in the distribution of mass 
amounts lost in the explosions; since the eccentricities for 
symmetric explosions is directly proportional to this SN 
mass loss, the double-peak signature is imprinted on the 
eccentricity distribution too. Once kicks are introduced 
(even of small magnitude) this structure is "washed-out" 
due to the randomizing effect of kicks. It turns out that 
this behavior has no effect on our statistical analysis that 
follows and therefore we do not discuss this in any further 
detail. 

3.2. Birth vs. Present Distributions 

We next compare the DNS eccentricity distributions 
(using models with Arzoumanian kick distributions) at 
two different times in their evolutionary history: at birth 
(the time when the second NS is formed) and at present 
(shown for three of our models in Fig.|3l2j), i.e., at 10 Gyr 
since the formation of the first primordial binaries in the 
simulations. Emission of gravitational radiation and re- 
sultant orbital inspiral, circularization, and mergers for 
the most short-lived DNS change the binaries' orbital ec- 
centricities as time progresses. We consider the present 
distributions a more realistic representation of the Galac- 
tic DNS population, since binaries must have experienced 
this kind of orbital evolution after their birth. In our sim- 
ulations, each primordial binary is assigned (randomly) 
a formation time, and is subsequently tracked at regular 
time intervals so that it is possible to know the binary 
properties binary at a particular time. 

It is evident that, for coalescing binaries (Fig.[3J), there 
is a noticeable qualitative change of the present distri- 
butions towards lower eccentricities. This is in agree- 
ment with suggestions m ade and results obtained by 
|Ch aurasia &: Bailesl l|2005|h On the other hand, for wide 
DNS (Fig.^J and for most binary evolution models, there 
is no clear evidence for significant evolution of the eccen- 
tricity distribution with time. We assert that this differ- 
ence is due to the fact that the ages, i.e., time elapsed 
from their birth to the present, of most DNS binaries 
are short compared to the typical timescale relevant to 
orbital evolution due to gravitational radiation. Clearly 
this is consistent with the fact that these systems have 
such properties that do not coalesce within a Hubble 
time. 

4. STATISTICAL ANALYSIS 

In what follows we present the im plementatio n of a 
Bayesian statistical analysis, following Finn ( 1994). This 
method allows us to compare the model predictions with 
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Fig. 3. — Eccentricity distribution functions of coalescing binaries 
for three of the models considered in our study, using Arzoumanian 
kicks. DNS populations are shown at birth (dotted lines) and at 
present (lOGyr, solid lines). 




eccentricity 

Fig. 4. — Eccentricity distribution functions for wide binaries. 
Line types as in Fig. 3. 



the small-number observed samples, and calculate the 
likelihood of each of the models we consider here. 

4.1. Method 

4.1.1. Observed Eccentricities 

A priori it is not clear that the measured DNS ec- 
centricities are indeed the true values, since measure- 
ments are associated with errors (no matter how small). 
Therefore, we assume that the measured eccentricities 
are distributed normally (following a Gaussian distribu- 
tion) about the true value of each binary's eccentricity: 



exp 

P(e\e,I) = 




Here, e corresponds to the observed eccentricity, e is the 
true eccentricity, and a e is the reported uncertainty in 
the measurement. / is a parameter that encompasses all 
the unspecified assumptions that are also inherent in the 
observations. Note that we do not really know e, the 
true value; we only know e and <r e , which are the ob- 
servations and errors of the eight DNS in the observed 
Galactic sample (Table 0. Of course, given the small 



errors associated with these measurements, the proba- 
bility that the true value is very much different than the 
measured value is very small (lies in the tail of the Gaus- 
sian distribution). P(e\e,I) will be used to describe the 
probability of individual eccentricities in our Bayesian 
analysis, which is described in the next section. 

4.1.2. Bayesian Statistics 

In what follows we describe how we can use the ob- 
served sample and a set of theoretical eccentricity distri- 
butions to calculate the conditional probability that the 
observed sample is drawn from a particular theoretical 
model. In Bayes' law of conditional probabilities: 

P(M\D,I) = P f}^ P{M,I\ (4) 

where P(M\D, I) is the posterior probability, i.e., the 
probability associated with a model M, given the data 
set D; P(D\M,I) is the likelihood of the data set, given 
a model and we will denote by A(D). P(M,I) is the 
prior probability of the model, and P{D\I) a normalizing 
factor for P(M\D, I). This normalizing factor P(D\I) 
is a constant so that J P(M\D,I) dM = 1, and it is 
unknown at this point in the analysis. 

Concerning the theoretical models considered here, we 
have no specific prior bias against one or another. In 
other words, all models M are considered equally likely 
a priori (in the absence of data), so P{M, I) must simply 
be a constant. Since this prior probability is a constant, 
it can be absorbed into the unknown normalizing factor 
P(D\I), so that 

P(M\D,I) = C-A(D), (5) 

where C = P(M, I)/P(D\I), C makes P{M\D,I) a 
properly normalized probability density, and D in our 
case is the set of observed eccentricities, which we will 
denote as e s . 

Since the observations and property measurements for 
each DNS in the sample are independent, the likelihood 
of the complete set of measured eccentricities is equal to 
the product of the individual measurement likelihoods: 

A(e s ) = Y[P(e l \M,I). (6) 

i 

P(e-i\M, I) can be expressed as a convolution integral of 
P(ei\e,I) over all possible e values, given a model M: 

P( ei \M,I)= [ P(ei\e,I)x P(e\MJ)de. (7) 

P{e-i\e, I) is just the Gaussian expression given in § 4.1.1 
for the distribution of measured values given a "true" 
value, and P(e\M, I) is the distribution of "true" ec- 
centricities, which we have generated from the Monte 
Carlo simulations of DNS binaries for each model. Con- 
sequently, for each model, we have at hand all the ingre- 
dients we need to calculate the likelihood of the data set 
A(e s ), using eqs. (4) and (5). The relevant integral is cal- 
culated numerically using our model distributions. Since 
A(e s ) is strictly proportional to the posterior probabil- 
ity of the model P(M\D, I) (see eq. [3]), we do not have 
to calculate the constant C, and instead we can evalu- 
ate the various models based on their likelihood values 
A(e s ). The results of this evaluation is described next. 
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TABLE 3 
Likelihood Values 



Model 


Close 


Wide 


A birth a x 10 43 


A Arz b x 10 4a 




Ai d x 10 4;i 


A 50 c x 10 4;J 


A Mrth a x 10^ 


A Arz b x 10 2X 




A u) d x 10^ 


A 5 o e x 10™ 


A 


5.00 


57.46 


0.00 


0.00 


14.30 


1.34 


2.84 


3.94 


6.40 


3.20 


C 


4.42 


32.99 


0.00 


0.00 


87.11 


1.01 


2.31 


4.85 


7.05 


2.99 


El 


12.76 


35.20 


0.00 


0.00 


0.00 


0.99 


2.05 


1.27 


3.47 


9.17 


Gl 


4.11 


77.71 


97.75 


0.00 


85.93 


1.07 


1.55 


4.66 


4.97 


2.46 


G3 


9.28 


88.56 


0.00 


0.00 


139.96 


1.28 


1.09 


5.83 


6.42 


2.96 


H2 


6.70 


34.93 


0.00 


0.00 


0.00 


1.90 


2.30 


4.63 


6.10 


2.43 


LI 


4.32 


93.40 


91.08 


0.00 


52.67 


5.62 


2.44 


4.27 


5.65 


3.76 


L2 


5.64 


6.16 


0.00 


0.00 


0.00 


1.13 


2.37 


4.27 


7.65 


2.45 


Ml 


5.67 


77.61 


0.00 


0.00 


0.00 


1.42 


2.53 


4.27 


10.68 


3.81 


M2 


5.39 


52.64 


0.00 


0.00 


58.38 


1.69 


2.39 


5.26 


6.71 


3.43 



a A values at birth for Arzoumanian kick distribution. 
b A values at present for Arzoumanian kicks. 
C A values at present for zero NS kicks. 

d A values at present for a a = lOkms" 1 Maxwcllian kick distri- 
bution. 

e A values at present for a a = 50kms _1 Maxwellian kick distri- 
bution. 



4.2. Statistical Results 

The full set of calculated likelihood values for various 
models are shown in Table |31 We note that A values 
are not normalized, so their absolute values have little 
meaning. However it is their comparative values that al- 
low us to address the questions of interest to us: given the 
observed sample of measured DNS eccentricities, (i) are 
the birth or present populations more likely? (ii) which 
of the ten binary evolution models considered here more 
likely? (iii) are "standard" kick magnitudes or small or 
zero kicks more likely? To address each of these ques- 
tions we examine the likelihood (or odds) ratios for the 
models relevant in each case and compare these ratios to 
unity. 



4.2.1. Birth 



Present Distributions 



To quantitatively compare our results for the present 
and birth distributions we present the model odds ra- 
tios A presen t/Abirth as a function of the binary evolution 
models (Pig-EJ. It is evident that for the majority of the 
models, the ratios exceed unity often by factors of more 
than 2 or as high as 20. We note that for wide DNS the 
ratio values "hover" within a factor of 2 from unity indi- 
cating that there is little or insignificant evolution from 
birth to present for these wide DNS binaries (as already 
noted in § 3.2). 

We conclude that the models for coalescing DNS pro- 
vide clear evidence that orbital evolution due to gravita- 
tional radiation is important. When allowed to evolve, 
the roughly uniform eccentricity distribution at birth 
shifts to favor lower eccentricities at present. Four of the 
five coalescing DNS inhabit the lower eccentricity part of 
the distributions, so the shift towards lower eccentricities 
usually raises the likelihood of present model distribu- 
tions over those at DNS birth. This is in full agreement 
with the conclusions reached by l|Chaurasia & Bailcs 
l2005l) based on a different analysis and arguments. On 
the other hand the lack of strong bias in favor of the 
present distributions from the wide DNS models is un- 
derstood, given their long orbital evolution timescales. 

Two exceptions to high odds ratios in the coalescing 
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Fig. 5. — Model odds ratios for the present vs. birth distri- 
butions for each of the ten binary evolution models. The squares 
correspond to close DNS, the triangles to wide DNS, and the circles 
to all systems combined. 



systems appear for models El and L2 (these, respectively, 
assume very low common-envelope efficiency and large 
specific angular momentum carried by mass lost from the 
binary in non-conservative mass transfer phases). For 
model El (see Fig- EJ , the eccentricity distribution al- 
ready favors low e at birth, such that when allowed to 
evolve, much of the DNS population possesses e < 0.05, 
a region where there are no observed DNS. We also see 
that the probability of seeing DNS at 0.2 < e < 0.3 ac- 
tually decreases, contributing to lowering the "present" 
likelihood. In model L2, the odds ratio close to unity is 
an indication of little evolution with time for the DNS 
formed under L2 assumptions combined with a negligible 
change in the PDF in the range 0.2 < e < 0.3. 

4.2.2. Preferred Binary Evolution Models 
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Fig. 6. — Model odds ratios for present distributions and for 
all models normalized to the standard model A: A preae nt,M values 



for all models, normalized to A 



close binaries, the triangles wide, 
wide binaries. 



Again, the squares mark 
and dots combine both close and 
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Fig. 7. — Model odds ratios, A ze ro/A y 4 rz for all models. The 
squares mark the close binaries; the triangles correspond to wide 
systems 



Guided by our conclusion that the present distributions 
are more likely, given the measurements, we compare the 
present likelihoods of the various binary evolution models 
to that of the standard model A (see Fig. EJ . 

It is interesting to note that for the majority of the 
models the ratios are well within a factor of 2 from unity. 
There is only one exception and this is model L2 with the 
assumed high specific angular momentum of mass lost in 
non-conservative mass transfer: its odds ratios for the co- 
alescing population compared to model A is lower than 
0.25. We conclude that, given the current size of the 
sample of measured DNS eccentricities, it is not possi- 
ble to distinguish between binary evolution models; the 
great majority of the ones considered here appear quan- 
titatively consistent with the eccentricity measurements, 
with the sole exception of L2. 

4.2.3. Zero and Low Kick Distributions 

Last we examine whether the origin of the low eccen- 
tricities is related to unexpectedly low or zero kicks im- 
parted to NS at birth. In Table |3 we list the likelihood 
values for all ten models and for models with zero or 
low kick-magnitude distributions at present. In Fig. [7| 
shows the model odds ratios for the zero kicks vs. the 
Arzoumanian kick distribution. 

In the case of wide DNS the ratios for most of the mod- 
els are within or very close to a factor of 2 from unity, 
with the exception of two models (Gl and G3). This 
could tentatively be used as an argument in favor of the 
hypothesis of very small NS kicks. However, it is most 
interesting to examine the case of coalescing binaries: the 
likelihood values are exactly equal to zero, for eight mod- 
els with zero kicks, for all ten models with kicks drawn 
from a a = lOkms -1 Maxwellian, and for four mod- 
els with kicks drawn from a a = 50kms _1 Maxwellian. 
As a result all these models are rendered highly unlikely. 
The reason for these values being identically zero can 
be traced back to the properties of the Hulse- Taylor bi- 



nary PSR 131913+16: it has an eccentricity of 0.617 and 
all the models above have distributions such that the 
probability for e > 0.5 is identically zero (see Fig. 0. 
Therefore, these models cannot allow for the existence 
of PSR B1913+16, and consequently their likelihood be- 
comes equal to zero. 

Based on the results for coalescing DNS, we conclude 
that models with vanishingly or moderately small kicks 
[a < 50kms ) are inconsistent with the current ob- 
served sample if our assumptions about the masses of 
exploding stars are correct. 

5. CONCLUSIONS AND DISCUSSION 

We address the question of the origin of the low 
eccentricities measured for observed DNS binaries 
(close and wide) that has attr acted considerable at- 
tention in the recent literature llvan den Heuveli r2004: 
[Ch aurasia fc Bailesl2 005L We use DNS population mod- 
els and associated predictions for eccentricity distribu- 
tions along with a Bayesian statistical analysis and quan- 
titative comparison with the observed sample. Assuming 
that all model DNS binaries can be detected as binary 
pulsars, we conclude that for the case of close (coalesc- 
ing) DNS binaries, models with vanishingly or moder- 
ately small kicks (a < 50kms _1 ) are inconsistent with 
the current observed sample (and specifically the Hulse- 
Taylor binary pulsar). We also examine the influence 
of orbital evolution due to gravitational radiation from 
the time of DNS birth to the present. We conclude that 
model distributions of DNS eccentricities at present are 
significantly more likely, given the observed sample. Our 
conclusion holds for the great majority of the models ex- 
amined, unless the birth properties are such that there 
is very little orbital evolution from birth to the present. 
These res ults are in agreement with the conclusions ob- 
tained bv lChaurasia fc Bailesl lj2005T) previously based on 
different methods of analysis. 

Here, we focus on just the distributions of eccentric- 
ities derived from simulations, since we are just inter- 
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Fig. 8. — The birth distributions of Tmr for models A, El, and 
M2, for two eccentricity ranges: 0-0.8 (dotted lines) and 0.8-1.0 
(solid lines). Arzoumanian kicks assumed. 



ested in the origin of the low values, and not to explicitly 
constrain binary evolution model assumptions. In forth- 
coming papers we will consider the much broader ques- 
tions of which binary models are most consistent with 
the complete set of observed DNS properties, such as 
orbital separations, masses, pulsar ages, in addition to 
eccentricities. Nevertheless we note that at a qualitative 
level, the range of orbital separations predicted by the 
StarTrack models for the DNS population is 1 — 6R© 
for the majority of coalescing DNS, in good agreement 
with the measured values; distributions of such DNS or- 
bital pro perties from StarTrack have been published in 

the past llBelczvnski et al.ll2002D. 

As also noted bv lChaurasIiTfc Bailesl l)2005l) . some frac- 
tion of the systems that have high eccentricities at birth 
actually end up merging very soon after (the rest of the 
high-eccentricity systems, evolve to lower eccentricities 
without necessarily merging completely by the present). 



As seen it is evident in Fig. [5] high eccentricity sys- 
tems at birth (0.8 < e < 1.0) tend to favor shorter 
merger times, more so than lower-e DNS. Highly eccen- 
tric DNS have accelerated orbital decay due to strong 
gravitational wave emission at periastron and due to 
their early mergers are not likely to be represented in 
the observed sample (as is the case for our Galaxy) . The 
importance of such short -lived systems has also been 
pointed out in the past iBelczvnski fc Kalogeral 120011: 
iGondek-Rosinska. Bulik. fe Belczvnskill2005|) . 

Current estimates of Galactic DNS merger rates (of 
interest to ground-based gravitational-wave interferom- 
eters like LIGO) are based on the DNS observed sam- 
ple and do not account for this selection effect against 
high-eccentricity DN S with short merger times (see 
iKalogera et alJ l)2004l) and references therein). Conse- 
quently these rate estimates could underestimate the true 
merger rate by possibly a significant factor. This is- 
sue has in fact point ed out in the past in a som ewhat 
different context (see iBelczvnski fc Kalogeral 1)200 lfl and 
llvanova et alJ l)2003|) ). To assess the magnitude of this 
upward correction factor that should in principle be ap- 
plied to such rate estimates, we calculate the fraction of 
DNS with merger times (at birth) shorter than 1 Myr 
and 10 Myr. We find that in most models these fractions 
are ~ 10 — 15% and ~ 20 — 30%, respectively; for a cou- 
ple of our extreme models (C and El) the fractions are 
~ 50% for a cut of 10 Myr. We conclude that although 
orbital evolution due to gravitational radiation affects 
the eccentricity distribution of the observed sample, the 
associated upwards correction factor to merger rate esti- 
mates is rather small (typically 10-40%). This factor is 
well within the current uncertainties of the DNS merger 
rate estimates for the Galaxy. 
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